home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / geom_lib / geomvals.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  9.4 KB  |  265 lines
  1. /******************************************************************************
  2. * GeomVals.c - Area, Volume, and counts on polygonal objects.              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March 1990.                          *
  5. ******************************************************************************/
  6.  
  7. #include <stdio.h>
  8. #include <ctype.h>
  9. #include <math.h>
  10. #include <string.h>
  11. #include "allocate.h"
  12. #include "convex.h"
  13. #include "geomat3d.h"
  14. #include "geomvals.h"
  15.  
  16. static RealType PolygonXYArea(IPVertexStruct *VHead);
  17. static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3);
  18.  
  19. /*****************************************************************************
  20. *  Routine to evaluate the Area of the given geom. object, in object unit.   *
  21. * Algorithm (for each polygon):                    V3         *
  22. * 1. Set Polygon Area to be zero.                   /\         *
  23. *    Make a copy of the original polygon             /      \          *
  24. *    and transform it to a XY parallel plane.           /        \V2         *
  25. *    Find the minimum Y value of the polygon           V4/         |         *
  26. *    in the XY plane                     \         |         *
  27. * 2. Let V(0) be the first vertex, V(n) the last one       \         |         *
  28. *    for i goes from 0 to n-1 add to Area the area of         \_______|       *
  29. *    below edge V(i), V(i+1):                     V0      V1      *
  30. *    PolygonArea += (V(i+1)x - V(i)x) * (V(i+1)y' - V(i)y') / 2             *
  31. *    where V(i)y' is V(i)y - MinY, where MinY is the polygon minimum Y value *
  32. * 3. Note that the result of step 2 is the area of the polygon itself.       *
  33. *    However it might be negative, so take the absolute result of step 2 and *
  34. *    add it to the global ObjectArea.                         *
  35. * Note step 2 is performed by another aux. routine below: PolygonXYArea.     *
  36. *****************************************************************************/
  37. double PolyObjectArea(IPObjectStruct *PObj)
  38. {
  39.     RealType
  40.     ObjectArea = 0.0;
  41.     MatrixType RotMat;
  42.     IPPolygonStruct *Pl;
  43.     IPVertexStruct *V, *VHead;
  44.  
  45.     if (!IP_IS_POLY_OBJ(PObj))
  46.     IritFatalError("Geometric property requested on non polygonal object");
  47.  
  48.     if (IP_IS_POLYLINE_OBJ(PObj)) {
  49.     return 0.0;
  50.     }
  51.  
  52.     Pl = PObj -> U.Pl;
  53.     while (Pl != NULL) {
  54.     /* Dont trans. original object. */
  55.     V = VHead = CopyVertexList(Pl -> PVertex);
  56.     /* Create the trans. matrix to transform the polygon to XY parallel  */
  57.     GenRotateMatrix(RotMat, Pl -> Plane);               /* plane. */
  58.     do {
  59.         MatMultVecby4by4(V -> Coord, V -> Coord, RotMat);
  60.  
  61.         V = V -> Pnext;
  62.     }
  63.     while (V != NULL && V != VHead);
  64.  
  65.     ObjectArea += PolygonXYArea(VHead);
  66.  
  67.     IPFreeVertexList(VHead);          /* Free the vertices list. */
  68.  
  69.     Pl = Pl -> Pnext;
  70.     }
  71.  
  72.     return ObjectArea;
  73. }
  74.  
  75. /*****************************************************************************
  76. *  Routine to evaluate the area of the given polygon projection on the XY    *
  77. * plane. Note the polygon does not have to be on a XY parallel plane, as     *
  78. * only its XY projection is considered (Z is ignored). Returns the area of   *
  79. * its XY parallel projection.                             *
  80. *  See GeomObjectArea above for algorithm:                     *
  81. *****************************************************************************/
  82. static RealType PolygonXYArea(IPVertexStruct *VHead)
  83. {
  84.     RealType PolygonArea = 0.0, MinY;
  85.     IPVertexStruct *V = VHead, *Vnext;
  86.  
  87.     MinY = V -> Coord[1];
  88.     V = V -> Pnext;
  89.     while (V != VHead && V != NULL /* Should not happen! */) {
  90.     if (MinY > V -> Coord[1])
  91.         MinY = V -> Coord[1];
  92.  
  93.     V = V -> Pnext;
  94.     }
  95.  
  96.     Vnext = V -> Pnext;
  97.     MinY *= 2.0;          /* Instead of subtracting twice each time. */
  98.     do {
  99.     /* Evaluate area below edge V-Vnext relative to Y level MinY. Note   */
  100.     /* it can come out negative, but thats o.k. as the sum of all these  */
  101.     /* quadraliterals should be exactly the area (up to correct sign).   */
  102.     PolygonArea += (Vnext -> Coord[0] - V -> Coord[0]) *
  103.             (Vnext -> Coord[1] + V -> Coord[1] - MinY) / 2.0;
  104.     V = Vnext;
  105.     Vnext = V -> Pnext;
  106.     }
  107.     while (V != VHead && V != NULL /* Should not happen! */);
  108.  
  109.     return ABS(PolygonArea);
  110. }
  111.  
  112. /*****************************************************************************
  113. *   Routine to evaluate the Volume of the given geom object, in object unit. *
  114. *   This routine has a side effect that all non-convex polygons will be      *
  115. * splitted to convex ones.                             *
  116. * Algorithm (for each polygon, and let ObjMinY be the minimum OBJECT Y):     *
  117. *                                V3         *
  118. * 1. Set Polygon Area to be zero.                   /\         *
  119. *    Let V(0) be the first vertex, V(n) the last.         /      \          *
  120. *    For i goes from 1 to n-1 form triangles           /        \V2         *
  121. *    by V(0), V(i), V(i+1). For each such           V4/         |         *
  122. *    triangle do:                     \         |         *
  123. *    1.1. Find the vertex (out of V(0), V(i), V(i+1))      \         |         *
  124. *         with the minimum Z - TriMinY.                 \_______|       *
  125. *    1.2. The volume below V(0), V(i), V(i+1) triangle,         V0      V1      *
  126. *      relative to ObjMinZ level, is the sum of:                 *
  127. *      1.2.1. volume of V'(0), V'(i), V'(i+1) - the                 *
  128. *         area of projection of V(0), V(i), V(i+1) on XY parallel     *
  129. *         plane, times (TriMinZ - ObjMinZ).                 *
  130. *      1.2.2. Assume V(0) is the one with the PolyMinZ. Let V"(i) and     *
  131. *         V"(i+1) be the projections of V(i) and V(i+1) on the plane  *
  132. *         Z = PolyZMin. The volume above 1.2.1. and below the polygon *
  133. *         (triangle!) will be: the area of quadraliteral V(i), V(i+1),*
  134. *         V"(i+1), V"(i), times distance of V(0) for quadraliteral    *
  135. *         plane divided by 3.                         *
  136. *    1.3. If Z component of polygon normal is negative add 1.2. result to    *
  137. *      ObjectVolume, else subtract it.                     *
  138. *****************************************************************************/
  139. double PolyObjectVolume(IPObjectStruct *PObj)
  140. {
  141.     int PlaneExists;
  142.     RealType ObjVolume = 0.0, ObjMinZ, TriMinZ, Area, PolygonVolume, Dist;
  143.     PointType Pt1;
  144.     PlaneType Plane;
  145.     IPPolygonStruct *Pl;
  146.     IPVertexStruct *V, *VHead, *Vnext, *Vtemp;
  147.  
  148.     if (!IP_IS_POLY_OBJ(PObj))
  149.     IritFatalError("Geometric property requested on non polygonal object");
  150.  
  151.     if (IP_IS_POLYLINE_OBJ(PObj)) {
  152.     return 0.0;
  153.     }
  154.  
  155.     ObjMinZ = INFINITY;     /* Find Object minimum Z value (used as min level). */
  156.     Pl = PObj -> U.Pl;
  157.     while (Pl != NULL) {
  158.     V = VHead = Pl -> PVertex;
  159.     do {
  160.         if (V -> Coord[2] < ObjMinZ)
  161.         ObjMinZ = V -> Coord[2];
  162.         V = V -> Pnext;
  163.     }
  164.     while (V != VHead && V != NULL);
  165.     Pl = Pl -> Pnext;
  166.     }
  167.  
  168.     ConvexPolyObject(PObj);           /* Make sure all polygons are convex. */
  169.     Pl = PObj -> U.Pl;
  170.     while (Pl != NULL) {
  171.     PolygonVolume = 0.0; /* Volume below poly relative to ObjMinZ level. */
  172.     /* We set VHead to be vertex with min Z: */
  173.     V = Vtemp = VHead = Pl -> PVertex;
  174.     do {
  175.         if (V -> Coord[2] < Vtemp -> Coord[2])
  176.         Vtemp = V;
  177.         V = V -> Pnext;
  178.     }
  179.     while (V != VHead && V != NULL);
  180.     VHead = Vtemp;       /* Now VHead is the one with lowest Z in polygon! */
  181.     TriMinZ = VHead -> Coord[2];         /* Save this Z for fast access. */
  182.  
  183.     V = VHead -> Pnext;
  184.     Vnext = V -> Pnext;
  185.     do {
  186.         /* VHead, V, Vnext form the triangle - find volume 1.2.1. above: */
  187.         Area = Polygon3VrtxXYArea(VHead -> Coord, V -> Coord,
  188.                                    Vnext -> Coord);
  189.         PolygonVolume += Area * (TriMinZ - ObjMinZ);
  190.  
  191.         /* VHead, V, Vnext form the triangle - find volume 1.2.2. above: */
  192.         Area = sqrt(SQR(V -> Coord[0] - Vnext -> Coord[0]) + /* XY dist. */
  193.             SQR(V -> Coord[1] - Vnext -> Coord[1])) *
  194.            ((V -> Coord[2] + Vnext -> Coord[2]) / 2.0 - TriMinZ);
  195.         PT_COPY(Pt1, V -> Coord);
  196.         Pt1[2] = TriMinZ;
  197.         if ((PlaneExists =
  198.          CGPlaneFrom3Points(Plane, V -> Coord,
  199.                     Vnext -> Coord, Pt1)) == 0) {
  200.         /* Try second pt projected to Z = TriMinZ plane as third pt. */
  201.         PT_COPY(Pt1, Vnext -> Coord);
  202.         Pt1[2] = TriMinZ;
  203.         PlaneExists =
  204.             CGPlaneFrom3Points(Plane, V -> Coord, Vnext -> Coord,
  205.                                     Pt1);
  206.         }
  207.         if (PlaneExists) {
  208.         Dist = CGDistPointPlane(VHead -> Coord, Plane);
  209.         PolygonVolume += Area * ABS(Dist) / 3.0;
  210.         }
  211.  
  212.         V = Vnext;
  213.         Vnext = V -> Pnext;
  214.     }
  215.     while (Vnext != VHead);
  216.  
  217.     if (Pl -> Plane[2] < 0.0)
  218.         ObjVolume += PolygonVolume;
  219.     else
  220.         ObjVolume -= PolygonVolume;
  221.  
  222.     Pl = Pl -> Pnext;
  223.     }
  224.  
  225.     return ObjVolume;
  226. }
  227.  
  228. /*****************************************************************************
  229. *  Routine to evaluate the area of the given triangle projected to the XY    *
  230. * plane, given as 3 Points. Only the X & Y components are considered.         *
  231. *  See PolyObjectArea above for algorithm:                     *
  232. *****************************************************************************/
  233. static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3)
  234. {
  235.     RealType PolygonArea = 0.0, MinY;
  236.  
  237.     MinY = MIN(Pt1[1], MIN(Pt2[1], Pt3[1])) * 2.0;
  238.  
  239.     PolygonArea += (Pt2[0] - Pt1[0]) * (Pt2[1] + Pt1[1] - MinY) / 2.0;
  240.     PolygonArea += (Pt3[0] - Pt2[0]) * (Pt3[1] + Pt2[1] - MinY) / 2.0;
  241.     PolygonArea += (Pt1[0] - Pt3[0]) * (Pt1[1] + Pt3[1] - MinY) / 2.0;
  242.  
  243.     return ABS(PolygonArea);
  244. }
  245.  
  246. /*****************************************************************************
  247. *  Routine to count number of polygons in given geometric object.         *
  248. *****************************************************************************/
  249. double PolyCountPolys(IPObjectStruct *PObj)
  250. {
  251.     int Count = 0;
  252.     IPPolygonStruct *Pl;
  253.  
  254.     if (!IP_IS_POLY_OBJ(PObj))
  255.     IritFatalError("Geometric property requested on non polygonal object");
  256.  
  257.     Pl = PObj -> U.Pl;
  258.     while (Pl != NULL) {
  259.     Count++;
  260.     Pl = Pl -> Pnext;
  261.     }
  262.  
  263.     return ((double) Count);
  264. }
  265.